MA/UC3M/10/96 



> 



a 



o 

CJ 



Closure of the Monte Carlo dynamical equations in the spherical 

Sherrington-Kirkpatrick model. 

L. L. Bonilla(*), F. G. Padilla(*), G. Parisi(**) and F. Ritort(*) 
(*) Departamento de Matemdticas, 
Universidad Carlos III, Butarque 15 



. Leganes 28911, Madrid (Spain) 

OO 

. E-Mail: bonilla@ing.uc3m.es 



E-Mail: padilla@dulcinea.uc3m. es 
E-Mail: ritort@dulcinea.uc3m.es 



fSJ ■ (**) Dipartimento di Fisica, 

o , 

. Universitd di Roma I "La Sapienza" 

0\ 



INFN Sezione di Roma I 
Piazzale Aldo Moro, Roma 00187 



C , E-Mail: parisi@vaxrom.romal.infn.it 



(February 6, 2008) 

Abstract 

We study the analytical solution of the Monte Carlo dynamics in the spherical 
Sherrington-Kirkpatrick model using the technique of the generating function. 
Explicit solutions for one-time observables (like the energy) and two-time 
observables (like the correlation and response function) are obtained. We show 
that the crucial quantity which governs the dynamics is the acceptance rate. 
At zero temperature, an adiabatic approximation reveals that the relaxational 
behavior of the model corresponds to that of a single harmonic oscillator with 



1 



an effective renormalized mass. 
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I. INTRODUCTION 



The off-equilibrium features of glassy systems constitute a very interesting field of re- 
search. A great deal of attention has been paid to the study of the dynamics of exactly 
solvable models of disordered systems where much relevant information can be gathered 
even in the mean-field case. Numerical simulations of short-ranged and long-ranged models 
reveal a similar nature of the off-equilibrum dynamics [|T]-|^. The outstanding phenomenon 
which characterizes the off-equilibrium regime is the presence of aging, which has been ex- 
perimentally observed in a large class of physical systems (for instance, in spin glasses |^ as 
well as in real glasses P). The distinctive feature of the aging phenomenon is that the re- 
sponse of the system to an external perturbation depends strongly on the time elapsed since 
the perturbation was applied. A phenomenological approach was proposed by Bouchaud 
where the slow dynamics originates from the presence of large number of traps in phase 
space with a very broad distribution of lifetimes. 

Much interest has been also devoted to the study of exactly solvable mean-field spin-glass 
models where information can be obtained from the solution of the dynamical equations. Of 
particular interest is the study of disordered models where replica symmetry is broken in the 
low-temperature phase. In this case, several results have been obtained in the study of p-spin 
models P|-pi][ as well as in the study of a particle in a random potential ]TT|-p!^]. From these 
studies a very interesting connection between mean-field glassy dynamical equations and 
the Mode Coupling Theory of glasses [Il5i|l6i has emerged. In this theoretical framework 
one derives a closed system of integrodifferential equations for the two-time correlation 
and response functions. Whereas all the information about the dynamics is contained in 
this system of equations, it is very difficult ascertain the long time evolution of one-time 
observables like the energy or the magnetization. Despite several advances in this direction 
(see for instance, the time series expansion carried out by Franz, Marinari and Parisi in 
case of the p = 3 spherical spin glass model [|1^) this is still an open problem. Recently 
a completely different approach to the dynamics of glassy systems has been proposed by 
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Coolen and Sherrington [|T^]. They introduce a closure assumption (the equipartitioning 
hypothesis) which yields evolution equations for the one-time observables. In comparison 
with the results of simulations they obtain fairly good results in the short time regime. 
However it is still unclear whether the equipartitioning approximation is applicable for long 
times and how to proceed in case it this is the case [IT^pn . 



On the issue of spin-glass dynamics, most of the existing works study Langevin equa- 
tions. Meanwhile Monte Carlo dynamics has received less attention notwithstanding its 



great importance in numerical simulations, p2|. In particular, given that a major part of 



the numerical work in spin-glasses uses the Monte Carlo algorithm, we think that it is im- 
portant to study the Monte Carlo dynamics in itself. Notice that Monte Carlo dynamics 
introduces the concept of acceptance rate which has no meaning for the Langevin dynamics. 
Interestingly, we shall see that the nature of the off-equilibrium regime strongly depends on 
the behavior of the acceptance rate. To shed some light on these issues we have consid- 
ered the spherical Sherrington-Kirkpatrick (SK) model While this model lacks replica 



symmetry breaking in the low temperature phase, it can be analytically solved. 

In a letter (hereafter referred as I) we introduced a method to solve the Monte Carlo 



dynamics of the spherical SK model [^]. We derived generalized dynamical equations for 
one-time quantities and a closed integrodifferential problem for their generating function. 
The same method can be successfully applied to models with Langevin dynamics. The 
purpose of this work is to present the method of the generating function, give a detailed 
account of the computations announced in I, and discuss some other topics such as the 
properties of two-time correlation and response functions and the long time behavior of 
dynamical quantities. 

The paper is organized as follows. In the second section we introduce the spherical 
SK model and explain how to solve its Langevin dynamics by the method of the generating 
function. While the results of this section are well known p^ , |25| , their derivation by means of 
the generating function serves as a clear introduction of this method. Section three contains 
the main result of this paper, i.e. the solution of the corresponding equations for the Monte 



Carlo dynamics. Section four analyzes the long-time behavior of dynamical quantities at 
finite and zero temperature and the conditions under which the results of the Langevin 
dynamics are recovered. The last section contains our conclusions while the Appendices are 
devoted to different technical matters. 



II. THE LANGEVIN APPROACH: A REMINDER 

The spherical Sherrington-Kirkpatrick spin glass model is defined by the Hamiltonian, 

^{^} = -J2 (1) 

i<j 

where the indices i,j run from 1 to (A^ is the number of sites) and the spins ai satisfy 
the spherical global constraint 

N 

E = N. (2) 

The interactions Jij are Gaussian distributed with zero mean and 1/A^ variance. The 
statics of this model reveals the existence of a thermodynamic second-order phase transition 
at T = 1 with a replica symmetric low-temperature phase. 



We will revisit the Langevin dynamics of the spherical SK model in order to introduce 
the technique of the generating function. The Langevin equation is 

_= + (3) 



with rii{t) a Gaussian white noise rii{t)rij{t') = 2T6ij6{t — t'), stands for average over 
the noise); fi{t) is a time dependent Lagrange multiplier which ensures that the spherical 
constraint eq.(0) is satisfied at all times and H is the Hamiltonian defined by eq.(|I|).We can 
rewrite the previous equation as a single mode problem, 

-^ = Xax + fi{t)ax + r]x{t) (4) 

Here a\ and rjx are the projections of the configuration {ai} and the original noise {rji} 
on the basis of eigenvectors which diagonalize the interaction matrix J^j. The {cta} satisfy 



the spherical constraint J2x=i c"a ~ "^^^ transformation matrix which diagonahzes the 
Jjj is an orthogonal matrix. Hence, the components rjx generate a white noise of the type, 



rjxit)'ny{t') = 2T6x,yS{t-t') . (5) 

The eigenvalues A of a random Gaussian symmetric matrix are distributed according to 
the Wigner semicircular law w{X) p7[1. 



^(A) = (6) 



A. Generating function for the one-time quantities 

Our purpose is to describe the time-evolution of certain one-time functions of the solu- 
tions of eq.@. We define the set of moments, 

/^. = ^E^^G^=^EA'^ (7) 

Notice that ho = 1 (spherical constraint) and hi = —2E where E is the energy. Using the 



result limt^t' Vxi^')'^xi^) = we get the equation 

dh 

= 2hk+i + + 2r < A'^ > (8) 

where 

«/(A)»= f f{\)w{\)d\ (9) 

and w{\) is given by eq.(^. In particular, eq.® for A; = gives the Lagrange multiplier as 
a function of the energy and temperature /i = 2E — T . In order to close these equations we 
define the generating function, 

1 1 CO k 

= ^^E^^e^^a, = ^E^'^^aW = T.T[hk{t). (10) 

A A;=0 

This function yields all the moments hk = ( ^ ll^fe'*'' ) q- 



We now want to formulate a problem for g{x,t) (equation, initial, boundary and sub- 
sidiary conditions) which has a unique solution. Then we will either be able to solve it (and 
thus determine the h'p.s exactly) or to infer the long-time behavior of the moments. By using 
Equations (|^, (||) and (|T0|), we find that g{x,t) satisfies the following differential equation 

dg{x,t) ^dg{x,t) 



dt dx 



+ 2/i(t)^(a;,t) + 2T<exp(xA) >, (11) 
M^) = -^(0,t)-T. (12) 



These equations have to be solved with the boundary condition g{0,t) = 1 (the spherical 
constraint) and an appropriate initial condition g{x,0) = go{x) (defined by the initial con- 
figuration {(Txit = 0)}). Besides taking the continuum hmit in the definition (p!0|), we see 
that 

g{x,t) = f^d\w{X)e^'g{X,t), (13) 

where the spectral transform g{X,t) > 0. The latter condition is kept by the dynamics if it 
holds initially. However its use is crucial to distinguish the physically meaningful stationary 
solutions from spurious ones. We can solve ([TT| ) by the method of characteristics assuming 
that /i(t) is given. The result is 

g{x, t) = go{x + 2t)e'Io ^(*')^*' + 2T /* dt' « e(^-+2{*^*'))A » ^2// ,it")dt" ^^^^ 

Jo 



Inserting (|Tj) in the spherical constraint g{0,t) = 1, it is possible to derive an integral 
equation which can be solved by means of the Laplace transform |2^ . 



It is easy to check that (^) is a stationary solution of eq.(|rT|) for T > Tc = 1. 



g-i{x) = - < — » (15) 



In this case the moments /i^ can be easily computed ||^ 



}f^ = -<^ ^ >, (16) 



and the Lagrange multiplier /x^'^ is given by 
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> 



For T < 1, ix'^'^ = —2, the stationary solution is given by 

1 



(17) 



(18) 



and the moments hi. are 



(19) 



Substituting eqs.(|T5| and [T^) in eq.(|TlD one can easily check that the right hand side of 
eq. fpUD is identically zero. 



B. The correlation and response function 

In order to find closed expressions for the correlation function we define the following set 
of two-times moments, 

1 ^ , ^ 1 



Note that in this notation the usual two-times correlation function is given by Co{t',t). 
The equation of motion for the Ck{t', t) reads, 



dt 



Ck+i{t',t) + fi{t)Ck{t',t) (21) 



where we have used the result r]\{t')a\{t) = for t' < t. Previous equation has to be 
solved with the initial condition Ck(t', t') = hk{t') where the hk{t) are the one-time moments 
previously obtained from the generating function g{x,t). 
We define the following generating function, 

1 1 r^k 

K{x,t',t) = ^Y.<t'){e-%a,{t) = -Y.^'-a,{t')a,{t) = Y.ljCu{t',t). (22) 

A fc=0 

The generating function K{x,t',t) yields the generalized two-times moments Ck{t',t) = 
( ^ ^dx^ ) Q ^^"^ satisfies the following homogeneous partial differential equation. 



dK(x,t',t) dK(x,t',t) . , 

with the initial condition K{x,t',t') = g{x,t') and g{x,t') is given in eq.(0). 

A similar method is applied to the response function. We define the set of two-times 
moments, 



G.{t,t) - - E(J - ^ EA ^ (24) 

where t' < t. In this notation the usual response function is given by Go (^'5 1). We construct 
the generating function 



Vix t' t)--^ ie^'^ ^^^('^) - 1 V .A.^^ 



00 k 

= T.l^G^it'^t)- (25) 
k=o 

The generating function r{x,t',t) yields the generalized two-times moments Gk{t',t) = 
( ^ ^dx'' ) 0' ^^^^ usual regularization of the response function at equal times the 
dynamical equation for the T{x,t',t) reads, 

"^^^^.'^''^^ = ^^^^'^''^^ + fiit)T{x, t', t) + 5{t - t') « exp(a;A) » (26) 
ot ox 

To solve this equation we need to impose the causality condition r(x, t',t) = if t < t'. 
The solution to equations (53) and (EBI) can be easily found. One gets the results, 



K{x, t', t) = g{x + t- t', t') el' 
T{x, t\ t) =< e("+*~*')^ > e/.' - t') (27) 

From these generating functions we can extract the usual two-time correlation and re- 
sponse function, Co(t',t) = K{0,t',t) and Goit',t) = T{0,t',t). At equilibrium we can 
substitute the solution (|l^) in (p6D obtaining, 

K{x, t\ t) = K^'^ix, t-t')= g^\x + t - t')e^''(*-*') 
r(a;, t\ t) = P«(a;, t - t') =< e("+*-*')^ > e^'^'^^'^'^Oit - t') (28) 



Both functions are translationally invariant (i.e. depend only on the difference of times) 
and the fluctuation dissipation theorem for the generating functions also holds, 

dK'^'^ix t — t') 

T^'^{x,t-t')=f3 '"' (29) 

The off-equilibrium behavior of the quantities Co(t', t), Goit', t) has been already studied 
in the hterature (see pl| , p5| ). 

III. THE MONTE CARLO APPROACH 

We consider the Monte Carlo (MC) dynamics with the Metropolis algorithm. The idea 
behind the MC approach is to postulate a dynamics in which a new configuration is proposed 
and accepted with a certain probability. The dynamics is ergodic at finite temperature and 
satisfies detailed balance. Following (I) we will consider a particularly simple motion which 
makes the dynamics exactly soluble: take the configuration {ci} at time t and perform a 
small random rotation from that configuration to a new one {rj} where 



and the are random numbers extracted from a Gaussian distribution p{r) of finite 
variance p, 

1 

p(ri) = ^=exp( V) • (31) 

Let us denote by AE" the change of energy Ai? = E{t} — E{a}. According to the 
Metropolis algorithm we accept the new configuration with probability 1 if /\E < and 
with probability exp{—(3AE) if AE > where f3 = ^ is the inverse of the temperature T. 

A. The probability distribution P{AE) 

As before in the Langevin formulation, we want to obtain one and two time moments. 
It is useful to work in the basis for which the interaction matrix Jij is diagonal. The energy 
then reads, 
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E{<y,} = -^Y.^<yl (32) 

where the ax are distributed according to eq.(|]). In this basis, the spin configuration still 
corresponds to a small random rotation, hence Cx ^ crx + rx/\/N where the new random 
numbers rx are extracted from the same Gaussian distribution eg. (|3T|) . 

The basic object we want to compute is the probability P{AE) of having a given variation 
AE of the energy E. This is a quantity which gives the average number of accepted changes. 
The variation AE* of the energy E in a Monte Carlo (MC) step is 



AE* = Xaxrx - — Y \' 



while the quantity ho = jr Sa^a is changed by the following amount. 



(33) 



^K^^^'T.n^J^^rl (34) 

The probability P{AE) of having a change of energy is 

P{AE) = [ S{AE - AE*)5{Ahl) l[{pirx)drx) (35) 

X 

where the last delta function in the integrand accounts for the spherical constraint and the 



variation AE* is given in eq-fpSf). 



Using the integral representation for the delta function 

1 r'^ ■ 

S(x) = — / e'^'^da (36) 

ZTT J-oo 

and substituting in (|35D we get 

P{AE) = j djjidri exp{iij,AE - 
^ ^a7a l>r-, n nAP^s\ /o7^ 

2iv ? ^rnf!y " 2 ? - -ir^) 

where 7a = /iA + 27]. After expanding the logarithm and retaining the first correction 
we get (after some manipulations) 
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with 

B, = h2-AE^- (/io = l); /^2 = ^E^'^a; (39) 

The equation for the energy is obtained by considering the average variation of energy 
in an MC step. In this case one MC step corresponds to N elementary moves. In the 
thermodynamic limit we can write the continuous equations, 



dE 
'dt 



/O foo 
{A)E P{AE)d{AE) + (AE) exp{~pAE) P{AE)d{AE) . (40) 
-oo JO 



We can check that equilibrium is a stationary solution of the Monte Carlo dynamics. 
Using standard static calculations p8|j21[] one can show that in equilibrium i?i = /i2 — 4i?^ = 



—2ET. In this case a straightforward computation shows that detailed balance is fulfilled. 
This means that, for a given value of AE, the first integral appearing in the r.h.s of eq.(^Dp 
cancels the contribution of the second integral in the r.h.s of eq. (^OD for the same value of 
AE. In other words, 

P{-AE) = exp{-l3AE)P{AE) (41) 



The equation for the energy reads, 

dE a{t) 



dt 2 

where the coefficients a{t) and b{t) are given by 



h2it) + h{t)E{t) (42) 



2 



a{t) = p2/5e 2 (/3Si+2i^) Erfip/SJ— - a) 

V Bi 

b{t) = -^{p'Erfia) + {AE - T)a{t)) (43) 

and Erf{x) is the complementary error function defined as Erf{x) = ^ dx exp[—x'^) 
and the parameter a is given by, 
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Note that the quantities E, Bi, and a depend on time. Also one can compute the 
acceptance rate as a function of time, which is the probabihty of accepting a certain change 
of the configuration, 

/O roo 
P{AE)diAE)+ exp{-PAE)P{AE)diAE) . (45) 
-oo Jo 

A straightforward computation shows. 

Ait) = + \e^iP^^^--) Erfipp^^ - a) (46) 



In order to obtain the time evolution of the acceptance rate and to solve the eq. (|42| ) of the 
time evolution of the energy we need to know the energy E and /i2 at time t. Unfortunately, 
one can see that the time evolution equation for h2{t) involves h^lt) and so on. This hierarchy 
of moments can be closed introducing a generating function |^ like has been done in the 
Langevin dynamics. This is the purpose of the next section. 



B. Generating function for the one-time quantities 

To close the equations of motion we consider the set of moments defined in eq.(^. The 
basic object to compute is the joint probability distribution P{Ahk, AE). This quantity can 
be written as, 

P{Ahk, AE) = f 6{Ahk - Ahl)S{AE - AE*)6{Aho) l[{p{r^)drx) (47) 

where the last delta function in the integrand accounts for the spherical constraint and the 
variation AE* is given in eq. (|33D while the variation Ah^. is given by, 

^K = 4^T.^'^^r, + j^Y.^''rl. (48) 

* A A 

Following the same technical steps as in the derivation of P{AE) we obtain the following 
result, 

P{Ahk, AE) = P{AE) P{Ahk\AE) (49) 
13 



where P{AE) is the probabihty distribution eq(|38D and P{Ahk\AE) is the conditional prob- 
abihty of Ah^ given AE. The final expression for the conditional probability is, 

1 , i^h + p\h- « ») + 2f (AE + p^E))^ . 
P{Ahk\AE) = exp( T^TTTri 

(50) 

with Cfc = h^k — hi; Bk = hk+i + 2Ehk; {ho = l]hi = —2E)] and the average <^ ... ^ has 
been previously defined in eq.(^. 

In order to obtain the dynamical evolution of the moments hk we have to compute 
its average variation in a MC step over the accepted changes of configuration. In the 
thermodynamic limit we can write the continuous equations, 

-± = Ahk= d{Ahk) Ahk 
at J-oo 

/O r oo 

d{AE)P{Ahk,AE) + d{AE) exp{-(3AE)P{Ahk,AE)) (51) 

The solution for a general integral of the previous type is shown in the Appendix A. The 
following result is obtained, 

dhkit) 



dt 



a{t)hk+i{t) + b{t)hk{t) + Ck{t) (52) 



where the time dependent quantities a{t) and b{t) are given in eq.(|^) and the coefficients 
Ck(t) are defined by 

Ckit) = {2E{t)a{t) - h{t)) < A'^ >= pM(t) < A'^ > . (53) 

Here A{t) is the acceptance rate defined in eq.(^). Note that the rate variation of the 
moment hk depends linearly on the moments hk and /ifc+i, but the coefficients a(t), h{t) and 
c(x,t) are nonlinear functions of hi = —\E and i?i = /i2 — h\ (second cumulant). It is thus 
reasonable to expect that the Monte Carlo dynamics is determined by the evolution of the 
first two moments. By means of the moments generating function g{x,t) of eq. ([lO|), we 
obtain from eq. (|5T|) . 
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^^^ = a{t)^ + b{t)g + c{x,t). (54) 

where the time dependent quantities a, b are functions of the two first moments E[t) and 
/i2(t) (whose relation to g is indicated below) defined in ( ^3]) and 

c{x, t) = {2E{t)a{t) - hit)) < e"^ >= pM(t) < e"^ > . (55) 

As in the case of Langevin dynamics, g{x,t) is a solution of eq. ( ^4] ) plus the following 
initial, boundary and subsidiary conditions: 

g{x,0)= go{x) 
giO,t)= 1 

(56) 

gU=o= -2E{t) 
1^ \x=0= h2{t) 

The second condition is the spherical one and the third and fourth define the first and 
second moments of the set hk- This linear partial differential equation equation can be readily 
solved with the method of characteristics. The general solution for a partial differential 
equation of the previous type ( ^4t) is shown in the Appendix B. 

C. The correlation and response function 

In order to find the dynamical equation for the set of correlation functions eq.(pOD we 
perform a similar computation as has been done for the moments h^. 

The elementary move eq.(|30|) at time t induces a change AE* and AC|; in the energy and 
the fc-moment of the correlation function Ck{t',t) (in what follows we take t' < t) defined in 
eq.®, 

1 ^ 1 



A^* = -^E^^A(t)rA(t) - ^^A(t) 



2 



AC.* = 47E^'^A(t')rA(t) (57) 
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The probability of having a change of the energy AE and the correlation ACjt is given 
by the joint probability, 

P(AE, ACfc) = / <5(AE - AE*)5(ACfc - AC^)5(A/io) \{{v{rx)dr^) 

Using the integral representation of the delta function, retaining only the terms of order 
and performing all the Gaussian integrals we find P{AE, ACk) = P{AE)P{AE\ACk) 
where, 

P{AE\ACk)= I ^=exp( 



where /i'a^ = h2k{t'y, Dk{t\t) = Ck+i{t\t) + 2E{t)Ck{t' ,t) and P{AE) is given by eq.®. 



To solve the equation of motion for the Ck{t',t) we write its average variation at time t 
over the accepted changes of configuration, 

dCk{t',t) 



dt 



/oo 
ACud{ACu 
-OO 



/O r oo 

d{AE) P(ACfc, AE) + / d{AE) exp{-(3AE) P{ACk, AE)) (58) 
- OO i/ 

Using the formulae of Appendix A we get the result, 

^.^|)c,„(*',0 + fc.(*',<) (59) 

where a(t),b(t) have been previously defined in eq. P3|) . 

Equation ( |5PD is solved with the initial condition Ck{t',t') = hk{t') (once the time evo- 
lution of the hk{t) has been obtained solving the preceding hierarchy corresponding to the 
set of moments hk). 

In order to close the previous hierarchy of equations we introduce the generating 
function K{x,t',t) of eq. (p2D which yields the generalized two-times moments Ck{t',t) = 
( ^ ^xfc* ) o' '^^^ K{x,t',t) satisfies the following partial differential equation, 

dK ait) dK b{t) 

-m = + 

16 



together with K{x^ t', t') = g{x, t'). The solution can be found by the method of the charac- 
teristics. A similar procedure can be used to obtain the generating function for the moments 
of the response function eg. (|2^) . In the Monte Carlo dynamics the equivalent quantity is 
given by the moments, 

= hm ^^A'=m,(t')^A(t) (61) 

A 

where A measures the intensity of an applied staggered field m^. This quantity measures 
the correlation between the spin configuration at time t and a small staggered magnetic 
field mx(t') applied at a previous time t'. Equation (|6l|) is the analogous of the correlation 
between the spin configuration at time t and the noise at time t' in the Langevin approach. 
The staggered magnetic field mx is an uncorrelated annealed random field taken from a 
Gaussian distribution of variance A, 



P(m.) = (27rA^)-texp(-^) (62) 



at each elementary move in the MC dynamics. The calculation of the response function 
proceeds in the following way: we compute the probability distribution for the variation of 
the Gk{t',t) when a small staggered magnetic field m\{t') is applied at an earlier time t' in 
an elementary move eq. (|30D at time t. We compute the joint probability distribution of the 
variation AGk{t',t) for an elementary move at time t, 

P{AE, AGk) = J S{AE - AE*)6{AGk - AGl)6{Aho) 

l[(p{rx)P{mx{t'))drxdmx{t')) (63) 

A 

where A/iq is given in eq.(|3^ and the variations of energy and G^ are given by, 
^E* = Y: AaA(t)rA(t) - ^ E ^A(t)' - mxit')rxit)) 

'A A 

^G^*fc = ^E^'^A(t')n(t) (64) 

Viv X 

Finally, the perform the limit A tending to zero for the intensity of the staggered magnetic 
field . Performing similar computations as in the correlation function case we obtain for the 
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average variation ^'^g* (in the region of times t' < t) the equation (|59[) substituting 
Ck by Gfc. Note that the difference between the hierarchy of equations associated to the 
correlation and the response function hes in the region of times t' > t where Gk{t',t) = 
while Ck{t',t) = Ck{t,t'). 

The partial differential equation associated to the generating function eq.(^) is in this 



case 



dT dT 

'dt ^ ar ^ ^'^^''^ ^ ^ exp(xA) > 5{t - t') (65) 
The partial differential equations (|60D and (|65|) can be also readily solved with the method 
of the characteristics (see Appendix B) as was done for in case of the one time quantities. 
Note that the equations for the generating functions g,K,T (eqs.(|^, |65|)) are formally 



the same as in the Langevin approach eqs.([TT|, |23|,|26D with the time-dependent parameters 
a(t), b(t), c{x, t) given by, 

aLANG{t) = 2- bLANcit) = 2{2E - T); CLANG{x,t) = 2T <^ exp{xX) (66) 

The main difference between the MC dynamics and the Langevin dynamics relies on the 
simpler time dependence of the coefficients a, b, c in the last case. This makes the large- 
time behavior of the Langevin dynamics exactly soluble (see P^ , p3[ ) while this is a very 
complicated task in the MC case. 



IV. ANALYSIS OF THE MONTE CARLO DYNAMICAL EQUATIONS 

In this section we proceed to solve the resulting dynamical equations for the MC dynam- 
ics. First we analyze the equilibrium dynamics showing that it coincides with the Langevin 
dynamics by an appropriate rescaling of time. Then we study the off-equilibrium behavior 
contained in the MC dynamics at finite and zero temperature. Since it is also our purpose 
to test the correctness of the solution of our equations we will compare the theory with real 
Monte Carlo numerical simulations. Moreover we will compare the resulting dynamics with 
that expected in the Langevin case. 
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A. Equilibrium Monte Carlo dynamics 



Dynamical equations for one-time and two-time quantities can be readily solved at equi- 
librium. Now the observables hk are independent of time and the two-time functions Ck(t', t), 
Gk(t', t) only depend on the differences of time t — t' . It is easy to check that the coefficients 
a{t),h{t),c{x,t) of eqs.p^E^) are time-independent and given by, 



p^l3{2E -T)Erf{a^'i) (67) 
d"i{x) = p2 < exp(xA) > Erf{a'"i) 

with 



. .lEL^ . Jil^ (68) 



This coefficients are the same as for the Langevin dynamics except for a rescaling of time 
t ^ t' = P l^Erf(a 1) ^^ Also ouc cau show (using eq.(^)) that the acceptance rate is given 



2 

by, 

A"' = Erf {a"') (69) 
Using the previous rescaling of time the equilibrium MC dynamical equations coincide 



with the Langevin ones [26 



B. Monte Carlo simulations 

In order to check our analytical results for the MC dynamics we have performed some 
MC simulations for finite sizes. To simulate enough large sizes we worked in the basis 
of eigenvectors cxx- In this way all the information about the quenched disorder is fully 
contained in the spectrum of eigenvalues A, which occupies much less memory than the full 
interaction Gaussian matrix Jij. The set eigenvalues A is chosen according to the semicircular 
law eq.(^. Typically we start from a random initial configuration (Xx = ±1 which fulfills the 

19 



spherical constraint. Then we perform a small rotation of this configuration {crx} {ctx + 
(eq.(pOD) where the are random numbers extracted from the Gaussian distribution 
eq.(|3TD. To keep invariant the spherical constraint we normalize the length of the vector 
{(Ta} in order to make it of length one. The resulting change of energy eq.(|32|) is computed 
and accepted with probability Min{l,exp{—i3AE)). We repeat this process for the new 
configuration and so on. A MC step corresponds to rotations. Because an elementary 
move involves a global change of the configuration {crx}, this algorithm is times slower 
than a usual MC algorithm with only local changes. We are able to simulate relatively large 
sizes in the range A^ = 500 — 2000 in a reasonable amount of computer time. 



In case of finite temperature it is relatively easy to solve MC equations in the large- 
time limit. All the information on the dynamics is contained in the time evolution of the 
coefficients a{t),b{t),c{x,t) which monotonically converge to their equilibrium values, then 
we expect the dynamics converges to the Langevin dynamics in the large-time limit except 
by a rescaling of the time 



At finite temperature and according to the magnitude of the parameter p we distinguish 
two different regimes depending if the acceptance rate is large or small. This two limits 
differ in the magnitude of the time scale above which Langevin behavior is recovered. 



In this case Erf {a) ~ 0(1) since a is small, i.e. the acceptance is always large. 
Note that in the particular limit p tending to zero the Langevin dynamics is recovered 
with a rescaling of time 



C. Finite-temperature dynamics 




(70) 



2 



1. The case p < 1 
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t - = (71) 
In figure 1 we show the decay of the energy for different temperature for a given value of 



/? = 0.1 by numerically solving the eg . (1631) . The theoretical prediction is compared with real 
Monte Carlo simulation results for = 500 spins. In the figure we also show the asymptotic 
large time-behavior in the Langevin case E{t) — E^'^ ~ (See ||2^) with the rescaling of 
time eq.([70|). 

Results for the correlation function Coityj, + 1) = -^'(O, tyj, + 1) are shown in figure 2 
for different values of = 1, 10, 100, 1000 at the temperature T = 0.4 with p = 0.1 obtained 
by numerically solving eg. (P4|) . Again we compare the solution with real MC simulations 
for N = 500 spins. In the figure we also show the asymptotic large-time behavior in the 

3 

Langevin case C{tw,tw + t) ^ (See P3[|). 



Note that in the dynamical regime shown in figure 2 there is no evidence of plateau or 
/9-relaxation process in the correlation function. This process is a signature of the existence 
of a short-time regime (FDT regime) where the system is locally in eguilibrium and the 
fluctuation-dissipation theorem applies (i.e. the Coit^^tw + t) ^ Co(t)). The reason why the 
/3-relaxation process is not seen in the figure 2 is that the values of shown there are too 
small. In the Langevin dynamics this regime usually appears for large values of which 
become much larger in the MC case since we have to rescale the time by the parameter 
eg.(|70D (the constant we have to rescale the time in figures 2 is approximately 87, i.e. 87 
MC steps correspond to one unit of time in the Langevin dynamics) . In this regime p <^ 1, 
because Erf {a) is finite, we expect the necessary to observe the plateau to scale like p~^ 
times the needed time in the Langevin case which reguires guite large values of tyj in the 
MC case. This off-eguilibrium regime observed in the case p <^ 1 where the acceptance rate 
is large is very similar to that observed in models with a continuous breaking of the replica 
symmetry, like the SK model [0] or finite-dimensional spin-glasses where no evidence 

of plateau is found in the region of values of explored. 
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2. The case p > 1 

In this case Erf (a) ~ since a is large, i.e. the acceptance is always very small. 
In this regime the rescaled time constant eq.(|70D is again small because Erf (a) is small (it 
decreases with p like pe~P^). For values of tyj ^ ^^^'^p ^ the system is in the short-time regime 
where a plateau is present and aging is present (i.e. the Co(t«,, + 1) strongly depends on 
tu))- This is clearly appreciated in figure 3 where we show the Co{tw,tw + t) for different 
values of t^, = 1, 10, 100, 1000 at T = 0.4, p = 5 and the numerical solution of eg. (P^) . In 
the regime < there is a plateau which increases with and eventually converges 

to the value qea = I — T for times of order This is trasient regime where the 

decay to the plateau does not satisfy time-translational invariance. No equivalent regime is 
found in the Langevin dynamics of the spherical Sherrington-Kirkpatrick model. 

The behavior shown in figure 3 is very similar to that found in systems with strong 
freezing in the low temperature regime where the acceptance rate is quite low. Generally 
these are systems with one step of replica symmetry breaking like the spherical or Ising 
p-spin model {p > 2) [|,^], the Potts Glass model with p-states {p > 4) [Q, the ROM 
(random ortoghonal model) and frustrated systems without quenched disorder like the 
Bernasconi model and other type of models |Q . All these systems are characterized by 
the presence of a dynamics with very low acceptance rate and the existence of a dynamical 
transition different of the statical one. While the statics of the aforementioned class of 
models is much different to that of the spherical SK model the dynamics shows interesting 
similarities. 



D. Zero-temperature dynamics 

In this subsection we are interested in extracting the large-time behavior of the dynamical 
equations at zero temperature. There is no reason in the zero-temperature case to get the 



same large-time behavior as in the Langevin case, the time-rescaling eq.([70D being ill defined. 
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Because in the zero-temperature regime the acceptance rate goes monotonically to zero (i.e. 
the parameter a diverges) the coefficients a(t),b(t),c{x,t) in eq.(^,p5|) have a non trivial 
large-time limit. Then we expect a dynamical regime much different from that predicted 
in the Langevin dynamics. In fact, we will see that there is no limiting case in which the 
Langevin dynamics is found. 

At T = the coefficients a{t),b(t),c{x,t) of eqs.(^3D, (|55|) become 

aT=o[t) - — 

&T=o(t) = -{Bt^ + '-^) (72) 
CT=o{x,t) = < e""^ > Erf (a) 
Also the equation (^) for the acceptance rate reads, 

Ait) = ^ (73) 

Note that the large time dynamics is governed by the time-dependent parameter a which 
diverges in the infinite-time limit. Unfortunately the dependence of the coefficients on a is 
strongly non-linear, which makes the mathematical treatment of the dynamical equations 
highly non trivial. However we can derive a nonlinear equation for the energy alone by 
means of the adiabatic approximation explained in the next subsection. Although a rigorous 
derivation of the adiabatic approximation is not known to us, the resulting equation yields 
an excelent approximation to the large-time behavior as given by the numerical solutions of 
the exact dynamical equations. 

We can analyze qualitatively how the system evolves at zero temperature. Suppose 
the system starts from a random initial configuration cxj = ±1 such that E{t = 0) = and 
i?i(t = 0) = 1. The energy monotonically decreases to the ground state energy E = — ii^s^^^ = 
— 1 while Bi decreases also to zero. In the large time limit a diverges and the acceptance 
rate goes to zero. There are two different regimes in the dynamics. The ffist one is an 
initial regime where a is small and the acceptance rate is nearly 1/2. This corresponds to 
a gaussian P{AE) (eq.(P^)) with width p\fB[ larger than the position of its center (p'^E). 
In this case, the changes of configuration which increase or decrease the energy have the 
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same probability. The energy decreases fast in this regime because the acceptance is large. 
The second regime appears when Bi is so small in order that a becomes large. In this 
case the acceptance is very small (it goes like and the dynamics is strongly slowed 

down. The system goes very slowly to the equilibrium. It is in this second regime where the 
adiabatic approximation developed in the next subsection applies. 

The adiabatic approximation shows that the large-time behavior in the Monte Carlo 
case is quite different form the Langevin case. This is not a surprise. Even for the simple 
harmonic oscillator, one can show that at zero temperature the dynamics is different in the 
MC and in the Langevin case [0. This is related to the non ergodic nature of the dynamics 



at zero-temperature in which only the changes which decrease the energy are accepted. 
This is clearly shown in figures 4, 5, 6, 7. In figure 4 and 5 we show the acceptation rate 
and energy as a function of time for three different values of p by numerically solving the 



off-equilibrium equations eq. (|B^) with the coefficients ([7^). 

We note two important results. First, as previously said, we observe that a is the 
parameter which governs the dynamics. More concretely figure 4 shows that the magnitude 
of a separates two different regimes. In the regime a <^ 1 the acceptance rate is 0.5 while 
it falls down rapidly to zero in the regime a ^ 1. Second, the large time behavior of the 
energy shown in figure 5 is strikingly different from that expected in the Langevin dynamics 
at zero temperature where E(t) + 1 ^ 0{l/t) ||25|. The energy decays like l/log{t) (plus 
corrections) as shown below. 

Figure 6 shows the CQ{tw,tyj + t) as a function of t for four different values of tw by 
numerically solving the off-equilibrium equations eq. ( |B4|) with the coefficients ([7^). Note 
that the large t behavior of the C{tw, t^+t) is strikingly different in the MC and the Langevin 
case since no trace of the power law decay is observed and the ^ scaley behaviour is 
lost. The correlation function strongly freezes (i.e. remains very close to one) for values of 
such that the acceptation rate at that time has already jumped to zero (see figure 4). This 
is reasonable because a very low number of accepted changes implies a very small change of 
the correlation function. 
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E. The adiabatic approximation 



We now discuss the adiabatic approximation which gives the correct asymptotic large 
time behavior of the energy and the acceptation rate at zero temperature. 

The equation of the energy eq.(^) can be written as a function only of the parameter 

a, 

This is an exact equation for the energy. Unfortunately we cannot solve it because we 
do not know the time evolution of a. Once a is known we also know the time evolution of 
the energy, hence also all the moments hk. In this sense, the parameter a, and consequently 
the acceptance rate, fully determine the dynamics. In the large-time regime, where a is 

_ 2 

very large, we can expand the error function Erf (a) ~ ^/^(l — ^) and we get the simple 
equation, 

= (75) 

E dt 2^a^ ^ ' 

We solve this equation using an adiabatic approximation which turns out to be the 
correct large-time solution as one can check numerically solving the equation ( |B3| ). We 
replace a by an effective parameter which depends on off-equilibrium quantities via a quasi- 
equilibrium relation. In our case we use the quasi-equilibrium relation eq. (^) with [3 = ^ = 
^j,*(E+i)) where n* is a renormalized parameter which has the physical meaning of an effective 
temperature. Note that /i* = 2 in the previous equation yields the equilibrium relation 
E = — in the low temperature phase. Substituting in (|75D and using Bi = —2fi*E{E+l) 
we obtain, in the large-time limit (i.e. ~ — 1) and neglecting higher orders in E + 1, 

In principle, /i* is an unknown parameter. In the simplest adiabatic approximation 
[I* = 2. This corresponds to equipartitioning in the surface of constant energy where the 
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equilibrium relation Bi = —4:E{E + 1) is fulfilled at all times. It is possible to show ||34 



that for /i* = 2 the equation ( |76D is the equation of the energy for a simple harmonic 
oscillator with Hamiltonian TC = ^muj'^x^ where muj'^ = A"*"^ = 2. Physically this means 
that the system relaxes as a simple harmonic oscillator with a value of mu;^ determined by 
the maximum eigenvalue of the spectrum. Physically a value of yU* < 2 means that the 
system relaxes as an effective single harmonic oscillator with effective mass mcu^ = larger 
than the single oscillator mass muj'^ = 2 . 

The large-time solution of eq.(^) can be easily worked out. One finds that the parameter 
a diverges like (log(t))5, the energy decays like E(t) = —l+0{j^^) and the acceptation rate 
goes like A(t) ^ j. All these quantities have non-trivial subdominant corrections. In figure 
7 we show the effective parameter /i* = — 2E{e+i) ^ function of time obtained by solving 
eq. (p3|) for different values of p. Note that the effective parameter /i* converges to a time- 
independent value in the asymptotic large time limit. By numerically solving the adiabatic 
equation ( |76D we have checked that the relaxation of the energy is in excellent agreement 
with the numerical solution of eq. ([B3D. This is a check of the adiabatic approximation which 
yields the values fi* ~ 0.087, 0.72, 1.65 for p = 0.1, 1, 5 respectively. Note that these quoted 
values of n* are smaller than 2. 



V. CONCLUSIONS 

In this work we have presented the analytical solution of the Monte Carlo dynamics of 
the spherical Sherrington-Kirkpatrick model. This is a very simple spin-glass model where 
full computations can be carried out. In order to analytically solve the MC dynamics we 
have used the technique of the generating function. The generating function allows us to 
close the dynamical equations in a similar way to how it has been done in other glassy models 
whithout quenched disorder. The present work is a step towards the use of this powerful 
technique in the spin glass context. As an example, we have solved the Langevin dynamics 
of the spherical SK model. In the context of the MC dynamics, the generating function has 
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been used to derive closed expressions for one-time quantities, like the generalized moments 
hk, and for two-time quantities such as the correlation or response functions. The formalism 
of the generating function may be used to compare the results obtained with both dynamics, 
MC and Langevin. 

By solving such a simple model, we lose the subtleties of having a replica broken low- 
temperature phase where the presence of a large number of metastable states makes dy- 
namics much more interesting |3^. Nevertheless our results for the MC dynamics make 
it clear how relevant for off-equilibrium dynamics is the acceptance rate, a quantity which 
is absent in Glauber or Langevin dynamics. Since a major part of the numerical work on 
glassy systems employs MC dynamics, its analytical study is of the outmost importance to 
be used as a guide when trying to extract conclusions from numerics. 

One of the most important results which emerge from this work concerns the importance 
of the acceptance rate in the off-equilibrium dynamics. We have seen that the main differ- 
ences between the Langevin and the MC dynamics appear whenever this rate is small and 
for times not too large. The Langevin dynamics is obtained in the limit p ^ with an 
effective rescaling of time t ^ t' = and the acceptance ratio is always 1 (cf. equations 
(^) and (^) ) . For finite p and finite temperatures we have seen that the Langevin regime is 
found for times t ~ p^^ when p is small and for times t ~ exp(p^) / p when p is large. In the 
intermediate regimes very many different behaviors are observed, specially for p > 1 where 
the acceptance is small. The dynamics observed in figure 3 reminds one a lot of what is 
observed in models with one step of breaking (see for instance ||3^j3^). Furthermore, at zero 



temperature, the MC dynamics yields a completely new behavior where the system rapidly 
freezes as soon as the acceptance goes to zero. No analogous regime is found in the case of 
Langevin dynamics where the transition from finite to zero temperature is smooth. 

We would like to comment three different directions which we think it would be inter- 
esting to explore. The first concerns the adiabatic approximation which is used to obtain 
closed equations for one-time moments. As it has been shown, this approximation suggests 
that the system behaves as a system of uncoupled harmonic oscillators with a renormalized 
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mass larger than the maximum mass corresponding to the largest eigenvalue of the Wigner 
spectrum. Although this suggestive approximation works extremely well, it would be inter- 
esting to understand it better. Future research should deal with a systematic derivation of 
the adiabatic approximation for one-time quantities (for instance, understanding the depen- 
dence of the effective parameter function of p), and its extension to the study of the 
long-time behavior of correlation functions. We think that the adiabatic approximation is a 
consequence of the relevance of entropic effects for the dynamics of the system. In the low 
acceptance regime, the system spends a long time searching configurations of lower energy 
in a very inefficient way. In this regime (not present when the dynamics is Langevin's) the 
adiabatic approximation works very well. More generally one would like to improve it in 
order to find a systematic way to close the dynamical equations without having to use the 



full generating function. The work done by Coolen, Sherrington and collaborators [T^ and 
that on the Backgammon model is a step in this direction. The second direction to 
explore is the study of models with one step of replica symmetry breaking (like the p-spin 
spherical spin-glass model |^) using the technique of the generating function. This is an 
interesting problem whose solution would allow us to obtain closed equations for one-time 
quantities like the energy and probably a set of higher moments. At present it is not clear 
whether this is possible or whether we can only obtain equations which relate the correlation 
and the response function Finally we want to mention that the type of differential 

equations studied here, (i.e. semilinear evolution equations like Eq. (Q) whose coefficients 
depend self-consistently on the first moments of the solution), appear often in the study of 
the dynamics of quite different models in statistical physics. Let us cite among others the 



problems of synchronization of populations of coupled oscillators |3^, arrays of Josephson 
junctions and plasmas or self-gravitating systems Often these problems can be 

described by non-linear causal equations where the time evolution of an observable at time 
t depends on all its previous history plus self-consistent relations for functions appearing in 
the equations. A general mathematical study of these types of equations would be welcome. 
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APPENDIX A: THE EQUATION OF MOTION 

In this section we write the general equation of motion associated to the general joint 
probability distribution 

P^AE,AO) = KeM-^=^)eM-^^) (Al) 

where c = e + f{AE — a) and a, 6, c, d, e, / are in general time-dependent parameters, AE 
stands for the energy change and AO the change of any observable (generalized moment, 
correlation or response function). In the MC dynamics a = —p"^ E, b"^ = p^Bi (with Bi 
given in eg. (p9|)) and c,d,e,f depending on the particular observable. The constant K = 
(47r6^(i^)~5 normalizes the probability distribution. 
The equation of motion for the observable O is, 

do f°° f°° 

— = / dx dyyP{x,y)w{x) (A2) 
where w{x) = Mm(l,e~^^) is the Boltzmann factor. Straightforward computations yield 



— = -{eErfia) + (e - fl3b') exp(-a/? + ^) Erf{^^ - a)) (A3) 



where a = — ^ and coincides with eq.(^4D with the previously quoted values of a,b. 
At zero temperature the equation of motion is. 



^ = l{eErf{a) - fJ^ exp(-«2)) (A4) 



with the same previous definition of the parameter a. 
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APPENDIX B: SOLUTION OF THE DIFFERENTIAL EQUATION 

In this section we present the general solution to the generating function g{x, t) in eq.i 
using the method of characteristics, 

^^ = a(t)|f + &(tk + c(x,t) (Bl) 

In order to solve this equation we make the following change of variables x ^ u{t) = 
X + /q a{t')dt' and eq.( |Bl| ) becomes 

'^^^ = h{t)g{u,t) +c{u,t) (B2) 
at 

which is a linear differential equation easily solvable. The final result is, 

g{x, t) = g,{x + a{t')dt')B{t) + B{t) |J dt>±±A0}^^ (B3) 

where B{t) = exp(/Q h{t')dt') and gQ{x) = g{x, t = 0) is the initial condition. Once the g{x, t) 
has been obtained one can also get the generating function of the two time quantities. For 
instance the K{x,t',t) associated to the correlation function in eq.(|60D is given by, 

K{x,t',t)=g{x + y'^ a{t")dt" ,t')[^f . (B4) 
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Figure Captions 

Fig. 1 Relaxation of the energy E{t) — E^'^ as a function of time for p — 0.1 and temperatures 
T — 0.2, 0.4, 0.6. The points are the MC simulations data {N — 500, one sample), the 
lines are the analytical solution of the MC dynamical equations and the dashed line 
is the large time behavior described by the Langevin dynamics rescaling the time as 
described in the text. 

Fig. 2 Co{t^, + 1) with p = 0.1, T = 0.4 for t^j = 1, 10, 100, 1000. The points are the MC 
simulations data {N = 500, one sample), the hues are the analytical solution of the 
MC dynamical equations and the dashed line is the large time behavior described by 
the Langevin dynamics rescaling the time as described in the text. 

Fig. 3 Co{ty,, ty,+t) with p = 5, T = 0.4 for ty, = 1, 10, 100. The points are the MC simulations 
data {N — 500, data averaged over 5 samples) and the hues are the analytical solution 
of the MC dynamical equations. 

Fig. 4 Acceptance rate A{t) at T = as a function of time for p = 0.1,1,5 obtained by 
solving the MC dynamical equations (lines) compared with MC simulations (points) 
(^N — 500, one sample). 

Fig. 5 Energy E{t) + latT = Oasa function of time for p = 0.1, 1, 5 obtained solving the 
MC dynamical equations (lines) compared with MC simulations (points) {N = 500, 
one sample). 

Fig. 6 Co{t^,ta, + t) at T = with p = 0.1 for = 1, 10, 100, 1000 (from bottom to top) 
obtained solving the MC dynamical i equations. 

Fig. 7 Ratio —Bi/{2E{E + 1)) as a function of time at T = for different values of p. The 
large-time limit yields the effective parameter p*. 
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c ( t „ + 1 ) 




